Logistička regresija

Predavanje VII

Autor

Doc. dr Nemanja Batrićević

Prije nego počnemo sa radom, potrebno je da učitamo prethnodno instalirane R pakete koji će nam biti potrebni za rad :

setwd("put/do/radnog/foldera/")

library(readxl) # paket za čitanje baze u .xlsx formatu
library(psych) # paket za opisivanje podataka (deskriptivna statistika)
library(car) # paket za rekodiranje
library(effects)
library(tibble)

dem <- read_xlsx("demlog.xlsx")

Učitavanje baze podataka demkat. Baza podataka sadrži informacije o državama:

Istraživačko pitanje:

Hipoteze?

Logistička regresija

Četvrti tip analize kojim ćemo se baviti na ovom predmetu odnosi se na statističke metode koje koristimo u situaciji u kojoj je zavisna varijable kategorička (dihotomna) a nezavisne varijable intervalne.

Ovaj tip analize pogodan je za sva istraživačka pitanja u kojima smo zainteresovani ispitati koja je vjerovatnoća da će se neki događaj desiti ili da će se neko ponašati na određeni način. U statističkom smislu, logistička regresija nam omogućava da damo odgovor da li prisustvo prediktora povećava (ili smanjuje) verovatnoću datog ishoda . Na primjer, kako ekonomske razvijenosti utiče na vjerovatnoću da država pripada kategoriji demokratskih država?

Ovo uputstvo pokriva slučaj kada je zavisna varijabla binarna (dihotomna) - to jest, đe može imati samo dvije vrijednosti, „0“ i „1“, koje predstavljaju ishode poput prolaska/neuspjeha, pobjede/poraza, ili zdravog/bolesnog. Slučajevi u kojima zavisna varijabla ima više od dvije kategorije ishoda mogu se analizirati multinomialnom logističkom regresijom ili, ako je više poređanih kategorija, ordinalnom logističkom regresijom. Međutim, ovim se analizama nećemo baviti na ovom kursu.

Prosta logistička regresija

Kreiraćemo ćemo model logističke regresije kako bismo predviđeli vjerovatnoću da je država konsolidovana demokratija na osnovu ekonomske razvijenosti. Funkcija glm() je skraćenica za generalizovani linearni model (GLM), koji pripada klasi modela koja uključuje logističku regresiju. Sintaksa funkcije glm() slična je sintaksi lm(), osim što moramo dodati argument argument family = binomial da bismo R-u sugerisali da sprovede logističku regresiju i da zavisna varijabla ima dvije kategorije.

S obzirom na to da je tip režima kategorička varijabla sa više kategorija, prvo je moramo rekodirati u dihotomnu.

table(dem$tip_rezima)

Authoritarian    Full Democ        Hybrid    Part Democ 
           52            24            39            52 
dem$konsdem <- recode(dem$tip_rezima, "'Full Democ' = 1; else = 0")
table(dem$konsdem)

  0   1 
143  24 
model1 <- glm(konsdem ~ oecd, family = "binomial", data = dem)

Logistička regresija suštinski procjenjuje vjerovatnoću da svaka država bude klasifikovan u određenu kategoriju - 0 (ostalo) ili 1 (konsolidovana demokratija). Vjerovatnoća se kreće od 0 do 100, pri čemu bi granica za klasifikovanje trebala biti logički povučena na 50%, jer u tom slučaju država ima identične šanse da bude u obje kategorije.

summary(model1)

Call:
glm(formula = konsdem ~ oecd, family = "binomial", data = dem)

Coefficients:
            Estimate Std. Error z value       Pr(>|z|)    
(Intercept)  -3.4012     0.5083  -6.692 0.000000000022 ***
oecdOECD      4.0431     0.6411   6.307 0.000000000285 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 129.522  on 152  degrees of freedom
Residual deviance:  72.704  on 151  degrees of freedom
  (14 observations deleted due to missingness)
AIC: 76.704

Number of Fisher Scoring iterations: 6

Slično linearnoj regresiji i logistički model možemo ocijeniti pomoću izvoda koristeći funkciju summary(). Imajte na umu da je format izvoda sličan onome koji smo viđeli u linearnoj regresiji; međutim, detalji o kvalitatu modela na dnu izvoda se razlikuju. Ovim ćemo se baviti kasnije, ali samo obratite pažnju na riječ devijacija. Ova devijacija (odstupanje) odnosi se na zbir kvadriranih odstupanja reziduala u linearnoj regresiji. Nulta devijacija predstavlja razliku između modela samo sa presjekom (što znači „bez prediktora“) i”zasićenog” modela (teorijski savršen model). Cilj je da devijacija modela (koja se označava kao preostala devijacija) bude što manja; manje vrijednosti ukazuju na bolji kvalitet modela. U tom smislu, nulti model pruža osnovu na osnovu koje se mogu upoređivati potencijalni model

Tumačenje koeficijenata

Tabela u nastavku prikazuje procjenu koeficijenta i ostale informacije koje proizlaze iz modela logističke regresije kako bi se predvidjela vjerojatnoća da je država konsolidovana demokratija pomoću članstva države u OECD. Ključni “problem” u slučaju logističke regresije jeste interpretiranje koeficijenata. Ono nije jednako intuitivno kao u slučaju linearne regresije. Koeficijenti u slučaju logističke regresije mogu uzeti tri zasebne forme: log-oddsrazmjer šansi (odds ratios) i procenti. Vodite računa da u originalnoj formi, procijenja vrijednosti koeficijenata iz logističke regresije karakterištu odnos između nezavisne i zavisne varijable na skali log-odds. Alternativne interpretacije zahtijevaju da transformišemo koeficijente.

Vidimo da koeficijent 4.04 sugeriše da članstvo u OECD povećava vjerovatnoću da je država konsolidovana demokratija.

Dalje možemo protumačiti koeficijent članstva OECD kao - razmjer šansi (odds ration). Rezultat ostaje isti, ali interpretacija zavisi od toga da li je koeficijent označen kao log-odds ili odds ration. Interpretacija je mnogo intuitivnija ukoliko koristimo razmjer šansi, pri čemu koeficijent 1.00 predstavlja referentnu tačku u kojoj je jednaka vjerovatnoća da osoba bude glasaš ili apstinent. Razmjer šansi dobijamo uzimanjem eksponenta od inicijalnog koeficijenta, koristeći funkciju exp().

exp(coef(model1))
(Intercept)    oecdOECD 
 0.03333333 57.00000000 

Na osnovu dobijenog koeficijenta (57.00) možemo reći da je vjerovatnoća da država članica OECD pripadati kategoriji konsolidovanih demokratija biti 57 puta veća nego u slučaju države koja nije članica OECD. U interpretaciji razmjera šansi, često se koriste probližne vrijednosti (npr. skoro šest puta veća šansa)

Mnogi aspekti izvoda modela slični su onima koje smo razmatrali u slučaju linearne regresije. Na primjer, možemo izmjeriti intervale povjerenja i tačnost procjene koeficijenta izračunavanjem standardnih grešaka. Na primjer, vrijednost p<0.000 sugeriše statistički značajnu vezu između ekonomskog razvoja (članstvo u OECD) i konsolidovanosti demokratije.

Kreiranje predikcija

Nakon što smo jednom procijenili koeficijente, na osnovu vrijednosti dobijenih u modelu možemo izračunati vjerovatnoću demokratske konsolidacije za karakteristike (varijable) koje smatramo relevantnim. Koristeći procjene koeficijenta iz našeg modela predviđamo da je vjerovatnoća za članicu OECDA da bude konsolidovana demokratija 65%. S obzirom da je ovo značajno veće od 50% možemo biti prilično sigurni da bi OECD država bila konsolidovana demokratija.

Ukoliko upredimo vjerovatnoće za obje grupe, vidimo da vjerovatnoća da je država konsolidovana demokratija opada sa 65% na 3% ukoliko država nije članica OECD-a.

predict(model1, data.frame(oecd = c("OECD", "Not OECD")), type = "response")
         1          2 
0.65517241 0.03225806 

Višestruka logistička regresija

Takođe, možemo proširiti ovaj model kako bi analizirali binarnu zavisnu varijablu koristeći više prediktora (nezavisnih varijabli). Možemo probati kreirati model koji predviđa vjerojatnoću da je država konsolidovana demokratija na osnovu članstva u OECD, nivoa ekonomske nejednakosti i nivoa pismenosti.

Rezultati su očekivani. Članstvo u OECD značajno doprinosti vjerovatnoći da je država konsolidovana demokratija. Takođe sa porastom nivoa pismenosti vjerovatnoća konsolidacije raste. Nakon što uzmemo ove dvije varijable u obzir, GINI koeficijent (ekonomska nejednakost) nema statistički značajan efekat. 

model2 <- glm(konsdem ~ oecd + pismen + gini, family = "binomial", data = dem)
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(model2)

Call:
glm(formula = konsdem ~ oecd + pismen + gini, family = "binomial", 
    data = dem)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -54.24775   25.56970  -2.122 0.033874 *  
oecdOECD      4.15749    1.15478   3.600 0.000318 ***
pismen        0.48206    0.23685   2.035 0.041818 *  
gini          0.10920    0.08514   1.283 0.199598    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 103.806  on 116  degrees of freedom
Residual deviance:  42.162  on 113  degrees of freedom
  (50 observations deleted due to missingness)
AIC: 50.162

Number of Fisher Scoring iterations: 10

Pretvorimo koeficijente u razmjer šansi (odds ratios) kako bi lakše interpretirali rezultate.

exp(coef(model2))
                      (Intercept)                          oecdOECD 
 0.000000000000000000000002757407 63.910741166387715850305539788678 
                           pismen                              gini 
 1.619400503452814943372572997760  1.115390789336615284810250159353 

Na osnovu rezultata višestruke logističke regresije sada možemo reći da članice OECD-a imaju preko 60 puta veće šanse da pripadaju grupi konsolidovanih demokratija. Sa druge strane, koeficijente koji se nalaze u intervalu od 1 do 2 možemo interpretirati kao procente. S obzirom na to, možemo reći da jedinica porasta pismenosti povećava šanse da država bude konsolidovana demokratija za oko 62%.

Kao i ranije, možemo veoma jednostavo kreirati prediktorski model kojim ćemo izračunati tačnu vjerovatnoću demokratske konsolidacije za državu određenih karakteristika. Na primjer, recimo da nas zanima vjerovatnoća demokratske konsolidacije za državu koja je nije članica OECD-a, ima osrednji nivo ekonomske nejednakosti i visok stepen pismenosti:

pred <- tibble(oecd = "Not OECD", gini = 50, pismen = 95)
predict(model2, pred, type = "response")
         1 
0.04778004 

Imajte u vidu da bi koristili funkciju predikt morate specifikovati vrijednost svih nezavisnih varijabli u modelu.

Država ovih karakteristika ima otprilike 5% šansi da bude konsolidovana demokratija. Kolika je vjerovatnoća ukoliko ostale karakteristike ostanu iste a GINI koeficijent poraste za 30 procenata?

pred2 <- tibble(oecd = "Not OECD", gini = 80, pismen = 95)
predict(model2, pred2, type = "response")
       1 
0.570516 

Sada je vjerovatnoća značajno veća, približno 57%.

Kvalitet modela

Do sada smo kreirali dva logistička regresiona modela i ispitali njihove koeficijente. Međutim, ostaje pitanje koliko su modeli zapravo korisni, odnosno koliko model dobro odgovara podacima. Osim toga, zanima nas i koliko su tačna predviđanja o tome u koju kategoriju određena država treba biti klasifikovana (konsolidovana demokratija ili ne)?

U linearnom regresionom modeli viđeli smo kako nas R-kvadrat može informisati o kvalitetu modela. Sada ćemo analizirati i nekoliko načina za procjenu kvaliteta logističkih modela.

Pseudo R-kvadrat

Za razliku od linearne regresije, ne postoji R-kvadrat statistika koja nam može jednostavno dati informaciju koji udio varijanse na zavisnoj varijabli objašnjavaju sve nezavisne varijable. Međutim, postoji niz Pseudo R-kvadrat statistika koji bi mogli biti od koristi. Najistaknutiji je McFadden’s R-kvadrat. Međutim, za razliku od R-kvadrata u linearnoj regresiji, logistički modeli rijetko postižu visoki McFadden R-kvadrat. U stvari, shodno riječima kreatora McFadden Pseudo R-kvadrat ≈ 0.40 predstavlja veoma dobro poklapanje sa podacima. McFadden’s Pseudo R-kvadrat vrijednosti možemo procijeniti koristeći paket pschl i funckije pR2():

list(model1 = pscl::pR2(model1)["McFadden"],
     model2 = pscl::pR2(model2)["McFadden"])
fitting null model for pseudo-r2
fitting null model for pseudo-r2
$model1
 McFadden 
0.4386741 

$model2
 McFadden 
0.5938427 

Vidimo da model1 ima nižu vrijednost (44%) od modela2 (59%). Ipak, vidimo da oba modela imaju zadovoljavajući kvalitet modela.

Validacije predikcija

Kada se razvijaju modeli za predviđanje, najvažnija statistika odnosi se na to koliko dobro model predviđanja vrijednosti zavisne varijable među obzervacijama izvan našeg uzorka. Prvo, trebamo koristiti procijenjene modele za predviđanje vrijednosti na našem skupu podataka. Kada koristite predviđanje, obavezno uključite argument type = response tako da predviđanje daje vjerovatnoću da će ispitanik/država biti klasifikovan u kategoriju “1”.

test.predicted.m1 <- predict(model1, newdata = dem, type = "response")
test.predicted.m2 <- predict(model2, newdata = dem, type = "response")

Sada možemo uporediti predikcije zavisne varijable sa mjerenim vrijednosti za svaki model i viđeti koji model najpreciznije klasifikuje rezultate. Možemo započeti korištenjem matrice (tablice) koja opisuje preciznost klasifikacije za svaki model. Svaki kvadrant tabele ima važno značenje. U ovom slučaju 0 i 1 u redovima predstavljaju jesu li birači glasali ili ne. FALSE i TRUE u kolonama predstavljaju da li je model predvidio da će birači glasati ili ne.

  • Tačni pozitivni (donji desni kvadrant): ovo su slučajevi u kojima smo predviđeli da će država biti konsolidovana demokratija i da zaista jeste.

  • Tačni negativni (gornji lijevi kvadrant): Predviđeli smo da države neće biti konsolidovana demokratija i one to nijesu.

  • Lažno pozitivni (gornji desni kvadrant): Predviđeli smo da će države biti konsolidovane demokratije, ali one zapravo nijesu.

  • Lažno negativi (donji lijevi kvadrant): Predviđeli smo da će države neće biti konsolidovane demokratije, ali one zapravo jesu.

Rezultati pokazuju da u slučaju model2 vidimo da je 14.5% država u bazi pripada grupi tačno pozitivnih, dok je 79.5% onih koje su tačno negativne. Dakle, ukupno 94% država je pravilno klasifikovano. Ostalih 6% obzeracija u našem uzorku pogrešno su klasifikovana na osnovu model2. Među njima, 4.3% lažno pozitivnih i 1.7% lažno negativnih.

list(
  model1 = table(dem$konsdem, test.predicted.m1 > 0.5) %>% prop.table() %>% round(3),
  model2 = table(dem$konsdem, test.predicted.m2 > 0.5) %>% prop.table() %>% round(3)
)
$model1
   
    FALSE  TRUE
  0 0.784 0.065
  1 0.026 0.124

$model2
   
    FALSE  TRUE
  0 0.795 0.043
  1 0.017 0.145

Kad uporedimo dva modela možemo viđeti da je procenat pravilno klasifikovanih obzervacija u model2 za 4-5% veći od pravilno klasifikovanih obzervacija na bazi model1.